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Abstract 

To study diseases with both vertical and horizontal transmission spread in a spatially heterogeneous environment, a multiple 
^1 species mathematical model for diseases involving both vertical and horizontal transmission is formulated. The model is based 
| on a network represented by a directed graph with n nodes, where nodes are locations, and links are outgoing or incoming move- 
ment flow between nodes. To study if the network can be invaded by the disease, an expression of the basic reproduction number 
( ' for the network is derived. The reproduction number is shown to be a function of both vertical and horizontal transmission 
^ . parameters. As an example, we determine bounds of the reproductive number for a Rift Valley fever (RVF) virus transmission 
' meta-population model. The expression of the exact value and bounds of the reproduction number can be computed for het- 
erogeneous networks with different weather conditions, different mosquito genera and species for different nodes. The roles of 
' the movement rates and disease parameters on the reproduction number are analyzed. The main contribution of our study is 
that we have derived an expression of the reproduction number for heterogeneous networks and have found the conditions when 
1 livestock movements play a role on the reproduction number. 
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1. Introduction 



Communicable diseases can be readily transmitted from one (region) to another (region) JH Q]. Transportation, travel, 

■ migration, and diffusion etc. have greatly changed the way that infectious diseases spread Jj]. Therefore, the temporal and 
^j. . spatial spread of infectious diseases ^ and the impact of population travel on it HI are considered to be very important flOl- 
\Q 1 The spatial spread of effective infection resulted from the introduction of infected agents was observed in different locations 

. at different times resulted in great economic losses, many animal and human cases, and even death. Noteworthy examples 

■ include smallpox in the New World in the 1500s yj] and plague in Europe in the 1300s QjJ, |4£] . More recently, HIV/AIDS in 1981, 
£\| . West Nile virus in North America in 1990s J^t], and SAS in Asia in 2003 |@] have shown that infections can spread over vast 
^ — I ■ regions, and even jump continents [Q]. 

The newly introduced diseases are a concern to both theoreticians and empiricists J3l. Some driving forces of the spatial 
spread have been well understood J3]. The spatial spread of communicable diseases is usually related to temporal-spatial move- 
ment of hosts in spatially heterogeneous environment [9]. On a local scale, it is through the contacts between susceptible and 

■ infected individuals [7]; On a relatively larger scale, it is through the contacts of individuals from different regions, connected by 
mobility. It is very complicated on a much larger scale such as transportation by airplanes |7j]. Mathematical modeling for com- 
municable disease studies has been playing an important role in policy decision making in many countries, such as the Unites 
States, Canada, Netherlands, and the United Kingdom [10]. Epidemiological modeling has the advantage of economy, clarity, 
and precision of mathematical formulation ifioll . Spatially structured models, such as the meta-population models or multiple- 
patch model have been widely used in epidemiological modeling to capture the effect of space 111 ill . Meta-population models 
describe an open sy stem, in which spatially discrete sub-populations connected by the movement of individuals between a set of 



patchy sites 112111311 . Multiple subpopulations are included in such models to provide an ideal framework under the assumptions 
that the host population can be decomposed into several patches |0] to study the dynamics of disease spread temporally and 
spatially 11411 . However, modeling the spatial propagation of communicable disease is complex A 1 511 - 

For a disease involving multiple species moving at different rates and/or different patterns, it is important to model the spatial 
dynamics accurately [7]. Geographic heterogeneities have been playing an important role in the spatial spread of communicable 
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diseases H. However, many models have been minimizing the geographic heterogeneities Q7[| . Very few exceptions exist 
and they were summarized by 0], One of the earliest exceptions considers spatial spread of influenza in the Soviet Union by 
partitioning the country into cities and incorporating migration and transportation. One approach adopted is to consider the 
mobility of individuals between discrete regions, considering that the transmission does not take place during traveling [15]. 
Then the situation is that of a directed network, with nodes representing locations and links representing the movement between 
the locations lfl5ll . The crucial role of keeping track of the mobility rates and movement patterns can be highlighted by foot and 
mouth outbreak in the UK in 2001 . The infected cattle were distributed to many farms before the movement ban was announced 
lfl6ll . Therefore, it is necessary to develop a transportation network to capture the spatial spread of foot and mouth disease 01 . 

Many communicable diseases are propagated by the two distinct mechanisms of both vertical and horizontal transmission 
Intl . Vertical transmission is the passing of infection from mother to a portion of offspring of the infected parentage Il7i[l8ll . 
Vertical transmission is often through eggs or seeds among insects or plants IU9I1 . A variety of diseases that transmit both 
vertically and horizontally have been discussed in |20ll . They include human diseases such as rubella, hepatitis B, Chagas' 
disease, and AIDS [20, 19]. It may be inferred that vertical transmission must be an important biological mechanism as seen 
from the prevalence of vertical transmission [20]. The role that vertical transmission can play in the transmission of infection 
is discussed elaborately in and the most important role is maintaining the infection 1221]. Several examples of vertical 
transmitted diseases are given and the mechanisms of vertical transmission are investigated in 12311 . The logical complement of 
vertical transmission is horizontal transmission. For animal and human diseases 11911 . horizontal transmission is often through 
direct or indirect contacts with infectious hosts or infectious vectors, such as mosquitoes, hosts, or other biting insects fl9Tl . 

In literature, theoretical analysis of the meta-population models is scarce |2fl. The advantage of analytic tools compared with 
purely numeric tools is that the analytic tools can provide prediction and analysis of the model behavior yl]. The reproduction 
number, Rq, is defined as the average number of new infected individuals produced by one infectious individual, in a population 
with only susceptibles 12411 . One can analyze the dynamics of an epidemic and determine the sensitivity of model parameters 
with the reproduction number. An expression of the invasion threshold, which is a function of the basic reproduction number, 
has been proposed by H25fl . The invasion threshold determines the minimum number of individuals traveling to have global 
infection [25]. A formula for the computation of the reproduction number for multi-species in multiple patches is derived in 
J3I]. It is arguably that the basic reproduction number Rq is the most important quantity in communicable disease modeling [24]. 
Theoretically, Rq plays an important role in the analysis, as well as insights from a infectious disease model 12411 . It is one of 
the most frequently used quantities to estimate the dynamics of emerging infectious diseases urgently during an outbreak and it 
provides insights when designing control strategies for established infections [24]. There are many methods of calculating the 



reproduction number. One of them is the next generation method developed by [26, 27] and popularized by 12811 . The method 



only includes the infected and asymptomatically infected states, which keeps the size of matrices relatively manageable 12911 . Ro 
is characterized by a process of producing new infections through transmission. Subsequent generations of infection growing 
denotes the growth factor per generation representing the potential growth. The growth factor is mathematical description of 
Ro \2A< 26ll . The original nonlinear ODE (Ordinary Differential Equations) system is linearized [24!] and the linear system of 
ODE is described by a matrix, which is usually called the Jacobian matrix [24]. The next generation matrix derived based on the 
Jacobian matrix relates the number of new cases in various generations 12411 . The next generation method provides the basis of 
defining and computing the reproduction number [24]. This method is especially widely used for compartmental models. 

To our knowledge, there is very little work on the reproduction number for meta-population with vertical transmission. The 
transmission dynamics of a parasite with both vertical and horizontal transmission has been modeled using only two differential 
equations ll30ll . In this simple situation, the reproduction number is the sum of vertical and horizontal transmission. However, it 
does not hold for a more complex simulation such as the model in [ 3 1 ] when the next generation matrix for vertical and horizontal 
transmission are not both scalars. It is hard to derive an expression form of Rq with both vertical and horizontal transmission for 
multiple species. A simple expression of the root of a third or higher order polynomial equation is hard to derive. The bounds of 
the reproduction number with respect to homogeneous populations for RVF virus transmission model considering both vertical 
and horizontal transmission homogeneous populations , have been presented by Ir3lll . 

We present the computation of the reproduction number and its bounds for compartmentalized models considering diseases 
with both vertical and horizontal transmission. This study considers meta-populations, in which the population is divided into 
discrete subpopulations and each of them is treated as well-mixing. Some authors call this type of models patch models. We 
assume that individuals move between different nodes and the disease can be transmitted within a node but not during traveling. 
In this paper, following yfl, a n— node S-E-I-R model with h species, among which g species transmit the disease both vertically 
and horizontally, and h — g species only transmit the disease horizontally have been formulated. It is a general model applied to 
network level. The model can be used to study the temporal-spatial spread of an infectious disease with both vertical transmission 
within the species and horizontal transmission between species. We derived the expression of the reproduction number which 
is a function of both vertical and horizontal transmission parameters. Finally, we compute the exact value and bounds of 
the reproduction number for RVF meta-population model as an example and analyze the factors that affect the reproduction 
number. The upper bound depends on both horizontal and vertical transmission, while the lower bound only includes horizontal 
transmission. 



2 



Summarizing the contribution of our work as follows. 



1 . Study the role of vertical and horizontal transmission on the reproduction number. 

2. Derive the expression of the exact value of the reproduction number for a general multiple-species meta-population model 
considering both vertical and horizontal transmission. 

3. Compute the bounds for Rift Valley fever meta-population model. 

4. Analyze the effect of movement and disease transmission parameters on the reproduction number. 

Our work can help compute the exact reproduction number and estimate its bounds easily, as well as evaluate the impact of 
vertical and horizontal transmission in disease spread. With our approach, the impact of movement on epidemic spread can be 
investigated. Mathematical derivation and numerical simulations have shown that the movement rates affect Rq for a system with 
heterogeneous disease parameters across the nodes. 

The paper is organized as follows. In Section [2] we describe the method of computing Rq with the next generation matrix 
approach and present the general model, and the derivation of the reproduction number. In Section [3] we apply the expression 
of the reproduction number to RVF meta-population model to compute Ro and its bounds and study the effect of movement, 
disease transmission parameters, and the size of network on the reproduction number. In Section |4] we summarize and discuss 
the simulation results and mathematical derivations. 



2. The reproduction number for diseases with both vertical and horizontal transmission 

2.1. The computation o/Rq based on the next generation matrix 

One of the frequently used methods is to compute the reproduction number as the spectral radius of the next generation 
matrix II271 13211 starting with the linearized system [[2411. For the computation, we only consider the compartments corresponding 



to infected individuals ll24ll . First, the original nonlinear ODE (Ordinary Differential Equations) system is decomposed into 
matrix & and "f , is the rate at which new infections appear in compartment i, and "V; - "Vr - with "fr denoting 
the transfer of individuals out of compartment i and representing the rate of transfer of individuals into compartment i 



112811 . Assume that the number of infected compartments is m, the m x m matrices F denoting transmission [24] and Vdenoting 
transition [24] are defined as: 

F = [— ^ ], V = [— ]. 

OXj OXj 

Matrices F and V are linearized matrices based on the original nonlinear matrices & and & and called Jacobian matrices, 
where Xj is the number or proportion of infected individuals in compartment j, j - 1, 2, 3, ... , m, and m being the total number 
of infected compartments, xq is disease free equilibrium vector. The (i, j) entry of F representing the appearance rate of infected 
individuals in compartment j produce infected individuals in compartment i 112811 . The (j, k) entry of V -1 represents the average 
time that an individual stays in compartment j, where i, j,k — 1,2, 3, ... , m 12811 . Hence, the (i, k) entry of F V -1 represents the 
expected number of infection in compartment i results from the infections in compartment k 112811 - The m x m matrix FV _1 is 
called the next generation matrix. Note that matrix F is nonnegative and V is a nonsingular M matrix. 



2.2. The meta-population model for diseases with both vertical and horizontal transmission 

We assume that the disease can be transmitted vertically and horizontally and only the species that transmit the disease both 
vertically and horizontally move with movement rates independent of disease status. The number of species that transmit the 
disease both vertically and horizontally and that of only that transmit the disease horizontally is g and h - g, respectively. The 
movement rates are assumed to be independent of disease status yfl. For some species, such as mosquitoes, it is not realistic to 



keep track of their origin and destinations 13311 . We only keep track of where the individual is at time t. The dynamics for species 



k (k — 1,2, • • -h) in patch i is given by the following system with 4hn + 2gn equations. 



dP ki 
dt 
dQki 



[(b ki (N ki -q ki Iki)-OkiPkMk) (1) 
[bmuhi - OkiQuWk) (2) 



dt 

dS ^ n n 

—ij- =dkiPki5(k) + b ki N k i5(k)-^p mki S ki I m ilN mi -d ki S k i+ ^ co kji S kj - ^ co k ijS ki (3) 



k=\ j=l,j*i J=hj*i 



3 



= 'y^mkiS kdmilN mi - S ki E ki - d ki E ki + ^ OJ k jiE k j - ^ (JJkijEu (4) 
k=l 7=1.7*' 7=1,7*' 

^ = OkiQkiS(k) + s ki E ki - y ki I ki + ^ oikjihj ~ ^ Witij/« (5) 

7=1.7*1 7=1.7*' 

= Jkihi - dkiRki + ^ ^kjiR-kj- ^ UkijR-ki (6) 
7=1.7*' 7=1.7*' 

Without loss of generality, we use the following indicator functions to distinguish the species that transit a disease both 
vertically and horizontally and those only transmit the disease only horizontally. 

MJ £ ><™ 

What we are considering for vertical transmission are oviparous animals. But the model can be modified easily and applied 
more widely. For the species that transmit the disease both vertically and horizontally, compartments P and Q representing 
uninfected eggs and infected eggs respectively are included for species with vertical transmission. The total number of eggs laid 
is buNu with b ki q ki I ki infected eggs and b ki Nu - b ki q ki I ki uninfected eggs. After development period, 6 ki P ki S(k) eggs develop into 
susceptible adults and 8 ki Q ki 6(k) eggs develop into infected adults. For the species that transmit the disease both vertically 
and horizontally, the number of newborn species k in location ; is b k ,N k i. The number of any species k being infected is 
Ylk=\PmkiS kiImi/N m i and he number of death from each compartment J is duJu (J=S, E, I, or R). After incubation period, 
s ki E ki transfer to infected compartment. Following infection period, juhi recover from the infection. The number of individuals 
in compartment J moving in and out of location i is £" =1 - +i oJkjihj and £" =1 - +i cj ki jj ki , respectively. 

2.3. The reproduction number for the network-based general model 

To determine the matrix F and V, the infected variables are ordered by compartment, species, and node, i.e. 



Qn,Qu, ■ ■ ■ Qu, QiuQn, ■ ■ ■ Qin, ■ ■ ■ Qa, Qa, ■ ■ ■ Q gn , 
E\\,E\2, ■ ■ • E\ n ,Ei\,Ei2, ■ ■ -Ein, ■ ■ ■ ,Eh\,Eh2, • • • E n „, 
hi, hi, ■ ■ -hn, hi, hi, ■ ■ -hn, ■ ■ ■ ,h\,hi, • •■hn- 



When multiple species are involved, the matrices become very large. To simplify the computation of the reproduction number, 
the matrices are decomposed into blocks to derive lower or upper diagonal matrices. The matrix F and V can be decomposed 
into blocks as follows, respectively. 



where 



Qgnxgn Ug n x2hn 



Olh, 



•nxgn 



w 2h , 



(©,-=! ®ki) Ognxlhn 

V H 



nxgn 



Qhnxhn ■^hnxhn 
Qfinxhn Qfinxhn 



V H = 



xhn 



Qhnxhn 



-®im=i^ * 



hnxhn 



U — [ Og„xhn ~ ©f = i(©; = i @ki) ®gnx(h-g)n ], 



W ■■ 



^-gnxn 
-g)nxn 
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Here the subscript of each block indicates the size of the block. The symbol is the matrix direct sum. i.e. A B — 



A 
B 



where A and B are matrices. Blocks M and X in matrix Vh can be decomposed into r blocks, respectively, i.e. M = © Mk 
and X = 0j =I Xk . Each block represents the k - th species. The matrix Mk is as follows. 



M k = 



dk\ + Ski + Yl)=i j -Mm 

-Uk\2 dfa + Sk2 + 2}=l,j5t2 ^klj 



-<*>kn2 



dkn + Skn + Z" = i Oiknj 



and 



X k = M k + tjjijki - Skd- 
1=1 

The block J[ in Fh can be written into a h x h block matrix and the (k, m) entry is. 

Si 



Matrix V is as follows. 



v- 



The next generation matrix FV 1 is as follows. 



FV' 1 = 

Since FV~ l is similar to the following matrix, 
/ 

. 0^(01,%) 1 







gnx2hn 



-^^(©^(©L ^)) t/v* 1 



-^^(©^(©L &)) ^ 



-^^(©^(©m I)) F H ltf 



/ 







F^" 1 - W(©f =1 (©" = i)^)^ 



Therefore, 



s I 

*o = P(FV- 1 ) = MFffV^ 1 - W(0(0 -WV H y ). 



(9) 



The complexity of computing Ro with Equation (O depends on the model. In the next section we apply Equation (0 to Rift 
Valley virus transmission model. 



3. The reproduction number for Rift Valley fever meta-population model 



The following RVF model is based on the model presented in [31]. Here we only take into account livestock movement 
between nodes. We computed the exact reproduction number and the bounds for the whole network. The parameters are shown 
in Table [TJ 

3.1. The compartmentalized meta-population model 

The general model in Section lZ2l is applied to study the dynamics of RVF virus transmission with r — 4, g — 1. The model 
is described in |31|. The parameter r, represents the recruitment rate of animals lf34ll . We only consider livestock movement 
here. In this model, we consider that transmission parameters are different for different nodes. The letter i in the subscript of the 
parameters or variables represents node i (i = 1, 2, ■ ■ • , n). 



dP 

-±=bu(Nu-qihi)-diiPu 
at 



(10) 



5 



-^^buquhi-OuQu (11) 
df 

—± = O u Pu ~ d u S u ~ PiuS uhilN 2 i (12) 
df 

d£i ■ 

— - = -dyEuNu/Ki +(3 2li S uhi/N 2i - e u E u (13) 
df 

-47 = OuQu - duhiNu/Ku + s u E u (14) 
df 

= 0h(Ph + Gii) - dxiS uNu/Ki (15) 

df 

-^=b 3i N 3i -8 3i P 3i (16) 
df 

— ^ = 6 3i P 3i - d 3i S 3i N 3i /K 3 -/32 3i S 3i l2i/N2i (17) 
df 

^ = -£ 3/ £ 3 ; - d 3i E 3i N 3i /K 3 + /3 23i S 3i hilN 2 i ( 1 8) 

df 

-^■=e 3i E 3i -d 3i I 3i N 3i /K 3 (19) 
df 
HAT,. 

— - = 03^3; - d 3i N 3i N 3i IK 3 (20) 
df 

AS " " 

—^=r i -PmS2ihilN\i-f} 32 iS2iI 3 ilN 3 i-d 2i S 2i + ^ ^j^zj- Z w o^2i 

(21) 

ip « « 

= PwSvIii/Nu +p 32i S 2i I 3i /N 3i - s 2i E 2i - d 2i E 2i + o)i E V ~ Z w ^ 2i (22) 

,, « « 

"jT = 2j Vjihj- 2_i ^ijhi - d 2i I 2i + s 2i E 2i - y 2i I 2i - n 2i I 2i (23) 

7=1, J*' j=hj±i 

dR " " 

— ^ = Z <4*2; - Z <4* 2i + r2l ' /2i ~ (24) 

dN n 11 n 

—jp = b 2i N 2i - d 2i N 2i - fi 2i I 2i + Y u2 ji S 2j- Z w y S « + X 

7=1.7*' 7=1 J^i 7=1.7*' 

n n n n n 

- Z ^ + Z 4^- Z Z 4^- Z (25) 

7=1.7*' 7=1.7*i 7=1.7*i 7=1>7*' 7=1.7*' 

^ = b 4i N 4i -/3 Ui S 4i Iu/Nu -/3 2 4ifS 4 ihi/N 2i - j3 34i S 4i I 3i /N 3i - d 4i S 4i (26) 
dE 4 i 

—r- = PuiS AihiINu +/3 24 ifS 4i I 2 i/N 2i +p 34i S 4i I 3 i/N 3i - e 4i E 4i - d 4i E 4i (27) 
dl 4i 

——- - e 4 iE 4i - y 4 [l 4 i - fi 4 il 4 i - d 4 il 4l (28) 
= jAihi - d 4i R 4i (29) 

df 

— = ^iA^i - AUiAi - ^4iA^4i (30) 



We apply the expression of the reproduction number in Equation (0 to RVF meta-population model to see the complexity of 
the computation for the model. 
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Table 1 : Parameters in the model omitting the node index 



Parameter Description 



Range or value Units 



Source 



Pn 
P21 
P23 
P32 
Pu 
Pi4 
P34 
1/72 

1/74 

1/di 
l/d 2 
l/d 3 
Vei 
I le 2 
l/e 3 
l/e 4 

q\ 

1/03 

Ki 

K 3 



contact rate: Aedes to livestock 
contact rate: livestock to Aedes 
contact rate: livestock to Culex 
contact rate: Culex to livestock 
contact rate: Aedes to humans 
contact rate: livestock to humans 
contact rate: Culex to humans 
recover rate in livestock 
recover period in humans 
longevity of Aedes mosquitoes 
longevity of livestock 
longevity of Culex mosquitoes 
incubation period in Aedes mosquitoes 
incubation period in livestock 
incubation period in Culex mosquitoes 
incubation period in humans 
mortality rate in livestock 

transovarial transmission rate in Aedes mosquitoes 
development period of Aedes mosquitoes 
development time of Culex mosquitoes 
carrying capacity of Aedes mosquitoes 
carrying capacity of Culex mosquitoes 
livestock movement rate from node i to node j 



(0.0021,0.2762) 1/day |35| 

(0.0021,0.2429) 1/day 1 35 ] 

(0.0000,0.3200) 1/day |35] 

(0.0000,0.096) 1/day El 

(0.00 1 , 0.002) 1 /day Assume 

(0.00004,0.00008) 1/day Assume 

(0.0005,0.001) 1/day Assume 

(2.5) 1/day S 

(4.7) 1/day |36] 
(3,30) days |35] 
(360,3600) days 1 35 ] 
(3,30) days |35] 

(4.8) days |35] 

(2.6) days |35] 
(4,8) days 1 35 ] 
(2,6) days \M 
(0.025,0.1) 1/day |35] 
(0,0.1) 1/day |35] 
(5,15) days 1 35 ] 
(5,15) days 1 35 ] 
10000 |35] 
10000 |35] 

(0 ,h) 



3.2. The reproduction number and its bounds for n (n ^ 2) locations 



where 



d_ 

dt 



d 
It 



Ey 
Iy 





'buquhi 







= 


bnqnhi 









b\ n q\nhn- 


1 x« 






lx4n 

PzuSuhi/Nzi 
P212S 12I22/N22 

P21nS lnhn/N2„ 

PmS 2 iIn/Nn + P321S 21/31/^31 

P\22$22l\2lN\2 + P322S22I32/N32 

P\2nS 2nhnlN\ n + /?32«S 2nhnlNi„ 
P23lS 3\hl/N2\ 
P232S 32^22/^22 

P23nS 3nhn/N2„ 

puiSnIu/Nn+fPuiSiihi/Nn +P341S41I31/N31 

^142^42^12/^12 + fP242S 42I22/N22 + P342S 42^32/^32 



fi\4nS4>J\n/Ni„ + fP24nS 4nhn /N 2n + P34nS 4n hn / N. 



in 



4nxl 
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Qr = 



BnQn 

QlnQln 



E-v = 



duE n N n IK\ + s n En 
duE n NvilK] + E l2 E i2 

d2iE 2 iN 2 i/K 2 + e n E 2 ] + E" =lj5tl oAE 2 j 
d 22 E 22 N 22 /K 2 + s 22 E 22 + 2" =lj¥2 ^pEij ■ 



2"=1,^2 W 2; £ 22 



d 2n E 2n N 2n /K 2 + E 2n E 2n + Yj"j^\Jtn ^jn^y - 2"=1,^„ ^IjEln 

d-i\E^\Ni\ I K-i + £31^31 

dyiEyiNyifKs + £32^32 

dlnElnNln/Ki + e^ n E^ n 
d\\E\\ + S41E41 
d\ 2 E\ 2 + £42 £42 

d4„E4„ + £4„ZS4„ 



-e n Qu +d n I n Nu/K i -euEu 
-#12612 + d\ 2 l\ 2 N\ 2 lK\ - s l2 E l2 

-6\nQ\n + d\ n I\„N\„l K\ - S]„E\ 2 



-s 2 \E 2 \ + d 2 ihiN 2] _IK 2 + y 2 \h\ + fiixhx + 
-s 22 E 22 + d 22 I 22 N 22 /K 2 + 722^22 + Hnhi + Ey=i,j#2 



7=1,7*1 w ji 



'2;- 



£7=1,7*2 



1,7*1 W l/21 



^2/22 



2 

7#i ^jn 



l,;¥« W «/2 



s 2 „E 2n + d 2n I 22 N 2 „/K 2 + y 2n I 2n + fi 2 „I 2 „ + £" =1 _ 
^31^31^31/^3 -£3i£'3i 

^32/32^32/^3 - £32^32 
dlnhnNln/Ki - Ei„E^„ 

dfthi - £41 £41 

^42^42 - £42^42 
d Aid An ~ SArJ^An 

Theorem 1. There exists a unique nonnegative solution of the total number of livestock in location i at disease free equilibrium 

Proof The total number of livestock in location i at disease free equilibrium is represented by [ A 7 ^, N% 2 ■ ■ • N% n J . To 
solve the total number of livestock in each node at disease free equilibrium, we need to solve the following equation system. 



,2 



-COX 



"12 

,2 



21 



^22 + £"=1,^2 
, ,2 



CO, 



J 2n 



2j 





\ Nil ' 




r\ 




N 22 




r 2 




- n;„ . 




- r„ . 



(31) 



Denote 



dn+XU^j 



"12 

,2 



da + Z"= 



, ,2 

1,7*2 u 2 j 



J 2n 



n2 
n-\ 



d2n + rf=[^ nj 



\ A^2 ■ • • A^in 1 * s a variable vector to be solved. We note that "W is a diagonal dominant matrix of its column entries 

IH], i.e. <Wn > YlLi.itj'Wij, for V/, where *Wy is the (i, f) entry of < W. By Theorem 1 in page 654 of Q, W is invertible. 
Hence, there exists a unique solution of Equation system (T3TT ). Moreover, by Theorem (O in the appendix, <W~ l is nonnegative. 
The unique solution is 

[N 21 ^22 ■•• ^] r = [^i ^ 2 °2 ••• Nl] T = < W- 1 [n r 2 ■■■ r n f 
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The partial derivative of and 'f' with respect to 



each variable at disease free equilibrium are as follows. 



Here 



U„ x &n 



V = 



®" =I #1; 0„x8« 



04nx4« MAnxAn 
^AnxAn ^AnxAn 



M 4n > 



®Anx 



— (®t= l (®;L l e */))4nx4« Xfaix4n 

U = [0 nx5n ®1 =1 (b u qu) 0„ x3 „], 

05nxn 

w = -(®uo u ) 

Qinxn 

Mi 

^ 3 

Ma 

^[5 ^6 ^[7 



(32) 
(33) 
(34) 
(35) 



5° S° S° S° 

M\ = ®" =l p 2 \i-+, M 2 = ®1=iPm— M-i = ®1 =l Pi2i—+, Ma = e" = i/?23i-7j, 

ly 2i ly 3i ly 2i 



s°. 

M 5 = ®" =l /3lAiJ^, M 6 = 0" =1 ^24i^, Ml = ®" =1 034ijj£- 

1/ 2/ 3/ 



The inverse of matrix V is computed as follows. 



y- 1 = 



n" _L 

/-hmrf 1 







where 



®t=i Af ifc 1 



The next generation matrix FV 1 is as follows. 



FV' 1 = 



By Equation (|9j, 



*o = p{F H V H x - W(0 ^)UV H X ) 



(36) 



(37) 



(38) 



Blocks M and X in Vh can be decomposed into 4 blocks represented by Mk and Xk, (k — 1 , 2, 3, 4), respectively. Matrices M, 
and Xj are related to Aedes mosquitoes, livestock, Culex mosquitoes, humans with k = 1,2,3,4, respectively. 



Mi 



dt\N n , „ 

Q ^12^12 







A', 



+ £12 












4^ 
^1 



+ Sin 
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M 2 = 



d 2 \ + £21 + 2"=2 ' 



-Wj 2 ^22 + «22 + Ej=l,j^2 w 2; 



,2 

,J «2 



2« 



<£>» + £2n + E"=l U 



n-i ,.,2 
"j 



X 2 = M 2 - e" =1 (y 2 / + m - e 2i ), 



(39) 



(40) 



M 3 



■ £31 








d 32 N 32 

IT +£32 



+ «3« 



x 3 = m 3 - e" =1 £3«, 



M 4 = 



t/41 + £41 

fi?42 + £42 














<^4« + e 4« 



X4-M4- e" =1 £4,-. 



Theorem 2. Based on the model presented in Section UJ] we have 



Proof By Equation d38l ). 
By (O, d, and © , 



p(F H V H l ) < «b «S p(F h Vh) + max(<7„). 



/?o = pCFV" 1 ) = pCFffVtf - W(0? =1 — )£/V#). 

PI/ 



where 



Y = 



-^(©"=l^)^H = 



3 «xn 03„ x3 « 



04nx4n 04 nx 4n 

Y Z 



$l =l {b u quWi 0„ x3 „ 

03«x« 03„x3n 



Since F^V^ and -W ©" =1 -g-UV^ both are nonnegative matrices by Theorem|6]in the appendix, 



We next want to show 



p(f h v h x )<p(fv- x ). 



p{FV- 1 ) p{F H V H l ) + maxfou). 



(41) 



We notice that X\ is a diagonal matrix, Xj 1 is also a diagonal matrix. Moreover, the eigenvalues of — W(©" =1 g-.WV^ are the 
diagonal entries of (®" =l (b u qiiW^ 1 . Hence, -W(©? =1 ^)UV H l = TOf 1 for some f. Here 



2) = 



04nx4« 04 nx 4« 
04«x4« Q<\nx4n 



and 



<2 = 



03nx3n 



3n 



0„x3« (®"^(b li quW l 
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From linear algebra, each column of P is an eigenvector of - W(©" =1 ^)UV^ . By direct calculation, 

P 

where 



3 \nx4n -C.AnxAn 



03«x 



^3«x3« 



-(^ = i(buque u ))X^M^ 0„ x3 „ 

03«xn 03„x3« 



X 



n x3« Inxn 
hnx3n 

3nxn 



Therefore 



F H V H l ~ Wi®U ^) uy H = FhV^P + D)P- 



We claim P 1 FhV h 1 P is a nonnegative matrix. 
First, 



p- y = 



<H~ X 
--C l JK l X" 1 



where 



^3n> 



03nxn ^3«x3n 



It is clear that "H , X , and -X VTW are all nonnegative matrices. Hence, P is a nonnegative matrix. So it is enough 
to show F H V H l P is a nonnegative matrix. In fact, 



F H V H X P : 



&X; l K =1 (®l 1 £ki)Wi 1 'H + m(X k [ )J ^(X-!)X 



where J{(X k )£ is a nonnegative matrix. Let M := J[X k (® 4 k=l (® n i=l s ki ))M k lr H + ^(X^'XJ. The left is to check that M is 
nonnegative. 



M 







MiX-^l^dM-' 











^-'(e^^iOM- 1 o WiX-Hei^dM- 1 o 

o ^[^'(©ti^OM,- 1 o o 

J{ 5 X- l (® n HX eu)M- 1 J{ 6 X 2 l (®U s 2i)M 2 l M 7 X- 1 (®l 1 e 3i )M- i 



BUibuquW- 1 0„ x3 „ 

03nxn ^3«x3n 
















^{4X2 




Jl(,X 2 







5 1 



s 1 



B" = i(^i^i' s ii)) x r lM r °»x3« 

03nxn 03„x3n 



The only possible negative entries are in (2, 1) and (4, 1) blocks. But the block in (2, l)-entry is 

MzX^l^dM^i^buq^X^ +J[ 2 X-\-® n i=l b u queu)X^ Mj" 1 = 

By assumption, X\ and M\ are both diagonal matrices. The last equality follows X^ l M^ 1 = Mj^Xj" 1 . Similarly, the block in 
(4, l)-entry is 



MsX^ieUsuW^ielfiuqidXi 1 +^ 5 ^r 1 (-©" =1 buqueidX^M^ = 



-i/vr>« 



1 /■ m'i 



Hence, FhV h 1 P is a nonnegative matrix. This proves the claim. Therefore 



By Theorem 2 in II 3 911 . 



pCFV" 1 ) = p{F H V H x - W{®U —)UV H l ) = p{P~ l F H V H x P + £>). 



p(FV- 1 ) < p{p- l F H V H l P)+p(D) = p(F ff V^) +p(D). 



(42) 



(43) 
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where 



p(D) =p(-W(®l 1 — WVj) = pdsUibuquW; 1 ) =p(©r =1 ?u) < max(q u ). 

&U i 

which is only related to transovarial transmission rate of Aedes mosquitoes. This finishes the proof. 



Next, we are going to give a bound of p(FnV H l ). By Equation d32b and d36b . 



FrrV, 



HV H 



J{X k \® A k=l {® n i=x e ki ))M- k l J{X k 



Then 



p{F H y H X ) = p{ftX k l {®U(®Uzki))M- k x ). 



(44) 



Theorem 3. Based on the model presented in Section UH] we have 



^jminixdpiX-'M- 1 ) ^ p(F H V H l ) ^ ^msxixdpQCfM?), 



where y = e " e2 ^ 12 ' /j2 "' + ^fe'fe 
wnerex, b u (f> u +eu) Mht+sx)' 



Proof Denote S := J[{®l =l (®'l x B ki )X k [ M k l ) 















M2X?(&l =1 e u )M^ ^X^i^sidM-i 







JUX 2 l (®l,s 2i )M- 1 











MsX^i^suW^ ^2 l ^l=i s ^ M 2 &iX?{®U s ^ M ? 
To compute the eigenvalue of S, we first calculate the characteristic polynomial of S as follows 



Si 

S 2 S 3 

S 4 

S 5 S 6 S 7 



H/-fl| 



,}/ -Si 
-s 2 ^/ -S 3 
-s 4 ^/ 

-S 5 -S 6 -By AI 

AI -Si 
= A n -S 2 AI -S 3 
-S 4 /I/ 



= A n \ 



I AS 2 l 
0/0 
0/ 



AI -Si 
-S 2 AI -S 3 
-S 4 /I/ 



= A n \ 

= A"\S 2 \ 
= A n \S 2 \ 
= A n \S 2 \ 
= A"\S 2 \ 



-/IS^S 3 



-S 2 




,i/ 



J 2 

-s 3 

,17 



-Si + A 2 S 2 l -^s;'s 3 



'2 

/I 



-s 4 



-Si + A 2 £- 1 -^S^Ss 



-s 4 

-S, + i 2 S 2 1 

-s 4 



/ 



-^.(SiS^ 1 -^ 2 s 2 - 1 s 4 1 +S 2 'S 3 ) 




-^(SiS^ - A 2 S 2 1 S 4 + Sj'S- ') Si - ^S" 1 



s 4 



^"|S 2 ||S 4 | I -^(SiS^ - ^s-'s- 1 + S 2 'S 3 ) 
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= ^ 2 "|S 2 ||S 4 | I A 2 S 2 l S 4 l - (SiS; 1 + S 2 l S } ) | 
= A 2n \S 2 \\S4\\S 2 l \\S 4 l \ | A 2 I n - (S4S2S1S; 1 + 
= I A 2 I n - (S4S2S1S, 1 + S4S3) I ■ 



4 

So matrix S has 2n zero eigenvalues. Moreover, the spectral radius of S is the square root of the spectral radius of 

■-I 



+ S4S3. By Equation (0D, we have 



p(*W) = y^fS^^Sl 1 + S 4 S 3 ) (45) 

= A /p(S 4 (S 2 Si+S3S4)S 4 I ) (46) 
= Vp(S 2 S! + S3S4), (47) 
Recall that y[ 2 , MuXi, Mi, JU, M3, X3 are all diagonal matrices. Then 

S 3 S 4 = (®'; =1 s 2i s 3i )^3X^M^Jl 4 X 2 l M 2 l . 

Denote 

BUBzfimfhli £2iS3iP32iP23i 

Afi 1 • (48) 

Oli(Ol + Eu) 03i(03i + £3i) 

Then 



mjnCtt)/^ 1 ^ *) < p(S 2 S, + S3S4) < maxOr^P^'Mj 1 ). 



Therefore, 



^rnin{ Xi ))p{X- l M- x ) < p^V" 1 ) < ^m^( Xi )p{X 2 i M 2 '). (49) 

This finishes the proof. 

Corollary: Suppose for V/, if the birth rate, incubation rate of mosquitoes, contacts rate from livestock to mosquitoes, and 
from mosquitoes to livestock are homogeneous for different nodes, i.e. 



bu = bu b^i-bi, £i, = £i, S2i=e%, ,£3; =£3, 
Pm-Pn, P2U = P21, Pm-Pn, Pm = Pi2- 

Then 



where 



p{F H V H l ) = ^XP(X 2 'M 2 1 ), (50) 



_ SlBlPnPll + £ 2 £3/?3 2 ^ 2 3 
X foi(Z>i+£l) ^3(^3 + ^3)' 



Proof: This directly follows Theorem[3] 



By (@2]>, (03]), and 1 

^XP(X 2 l M 2 l ) < pCFV -1 ) «S ^XP(X 2 1 M 2 X ) + max(q u ) (52) 

Theorem 4. Besides the assumptions in Corollary 1, we further assume thatforV i, the death rate, incubation rate, and mortality 
rate of livestock are homogeneous for all nodes, i.e. <f 2 ; = c/ 2 , £ 2 ,- = £ 2 , p.2i — P2- Then 



(d 2 + s 2 )(d 2 + J2+ y" 2 ) "V (d 2 + s 2 )(d2 + 72+^2) 



V < Ro < J-Tj \7T~ \ + (53) 
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Proof Denote by A, (A) the z-th eigenvalue of matrix A. Note that X 2 and M2 commute, then 

A i (M 2 X 2 ) = Ai{M 2 )Ai(X2). 

By Gershgorin circle theorem, 

min \Ai(M 2 )\ - d 2 + s 2 

i 

min \Aj(X 2 )\ = d 2 + y 2 +p 2 



Furthermore, 



Therefore, 



By Equation 



min \Ai(M 2 X 2 )\ = (d 2 + s 2 )(d 2 +y 2 + p 2 ) 



p(x 2 1 m 2 1 ) = 



1 



1 



min ; \Ai(M 2 X 2 )\ (d 2 + s 2 )(d 2 + y 2 + p 2 ) 



R H _ I £2 r Sij3uP 2 \ + £3^32^23 1 

y(b 2 + s 2 )(b 2 + y 2 +ij. 2 )h 1 (b 1 +e 1 ) b 3 (b 3 + e 3 V 



X 



(d 2 + s 2 )(d 2 + y 2 + ju 2 ) 



X 



(d 2 + s 2 )(d 2 + y 2 + p. 2 ) 



+ max(gi;). 



(54) 
(55) 
(56) 

(57) 
(58) 

(59) 
(60) 



If the transovarial transmission rate of Aedes mosquitoes are homogeneous, i.e. for V/, qn = q\, then the bounds of the 
reproduction number is the same as the bounds for local transmission. 



Example 1. Compute the bounds ofRa ifX 2 M 2 is positive or reducible. 

Since X 2 , M 2 are both diagonal matrices, by Theorem ©, M 2 1 and X 2 1 are both nonnegative matrices. But they are not necessary 
to be positive or irreducible. If we further assume X 2 l M 2 1 is a positive or nonnegative irreducible matrix, then we can compute 
the bounds of Rq for heterogeneous cases as follows. 
According to Theorem (f5J), 



1 



max,- (d 2i + s 2i )(d 2i + y 2i + p 2i ) 



^p{X 2 x M 2 l ) sC 



1 



min/ (d 2i + e 2i ){d 2i + y 2i + p 2i ) 



Therefore, 



min ; Or,) 



< P(Ro) < 



max,- (xd 



min,- (d 2i + s 2i )(d 2i + y 2i + p 2i ) + 



max^;) 



(61) 



(62) 



^ max,- (d 2i + e 2i )(d 2i + y 2i + p 2i ) 
Example 2. : For n — 2, we want to show the bounds ofRo depend on the movement if diseases parameter are heterogeneous. 
Matrix M-„ X„ Mr 1 , and XT 1 (i = 1, 2, 3, 4) are as follows. 



Mi = 



fli 
a 2 



, M 2 = 



Mr 1 = 



«3 



'12 



i 
- L 

"2 



Cl4 



, M 3 = 



as 
«6 



, M 4 = 



dn + E4 
d\ + £4 



(63) 



, Mr 1 = 





j_ 







M^ 1 = 





a 6 . 


, M4 1 = 








1 
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X, = 



a-; 
a 8 



, *2 = 



flg 
— fc>?. 



- W 21 
«10 



, *3 = 



an 
an 



, X 4 = 



d 4 
d A 



(64) 



^ o 

-i 

"8 



gig 



o 
i 



x: 



z 



Matrix £2 is defined as follows. 



1 - w 12 w 2I 

W 12 1 ~~ W 21 



where 



•fill 



«2 



«4 = fl?22 + £22 + &>2j 



+ £12 A3 = 021 + £21 + W 12 , 

, few?, , 

fl 5 = — jr^ + e 3 ] fl 6 = — £p + £32, 



«7 = 



ft 

d 12 JV° 



A 



OlO = d 2 2 + 722 +H22 + o> 2 , flu = 



flg = fl* 2 l +721 +m +a) 12' 



J 2\ 

Matrix X 2 l M 2 1 is a matrix related to movement rate. 



X?M? 



an = 



K 3 ■ 



The spectral radius of X 2 1 M 2 1 is as follows 



(fl3«4 - OJ^WjjXagaio - (O^OJ 2 ,,) 

fl3fl9 + fl4fli() + ^ 2 n u \\ + ^ 



2 2 2 2 

03^9 + OJyJjJ^x fl3W2[ + fli()W 9 [ 

2 2 2 2 

fl4Wj2 + flgW^ £Z4£Z 10 + W^^j 



(65) 



where : 



5 = (0309 + Ljf 2 a> 2l ) 2 + 4wj 2 ^2i( fl 3 + a io)( a 4 + ^9) + («4flio + ^^^l) 2 — 2(01309 + a>\ 2 (jJ 2i )(a 4 a\Q + wfjW^) (67) 



2(0304 — «,,<A>2,)(a9aio — arLwi.) 



2 x2 



(66) 



By calculation, the spectral radius of X 2 X M 2 1 is related to movement. Therefore, Rq is related to movement. 
3.3. Simulation results 

First, we construct a 100 node network with homogeneous parameters across all nodes except movement rates. The disease 
parameters for each of the 100 runs are uniformly distributed within their respective ranges [40], but they are kept homogeneous 
across all nodes. We assume the movement rate between any two nodes is uniformly distributed between and - in each run. 
The reproduction number Rq is computed according to Equation ( 1381 ). The upper bound and lower bound of Rq are computed 
according to the expression in d41l . The lower bound of Rq (denoted by Rq) vs. Rq in each run are shown in Figure [T(a)| and 
the upper bounds of Rq (denoted by R^) vs. exact Rq in each run are shown in Figure [T(b)| As is shown in the figure, the upper 
bounds is slightly greater than Rq and the lower bound is slightly smaller than Rq. Through extensive numerical simulations, we 
have observed that the closer the parameters between nodes, the tighter the bounds are. 

Second, we construct three networks with 3, 4, and 100 nodes, respectively. We run simulations with all parameters in Table 
Q] except livestock movement rate a> 2 j held constant during 100 runs for each network. The reproduction number Rq does not 
change during 100 runs for each network. Moreover, the values and bounds of Rq obtained by numerical simulations are the 
same for networks with 3, 4, and 100 nodes. Therefore, Rq does not depend on livestock movement or the number of nodes in a 
network when the parameters related to RVF are homogeneous across all the nodes. 

To see the impact of movement rates on the reproduction number, we run simulations for four node networks and explore the 
impact of various parameters in different scenarios listed in Table|2] In scenario 1, the contact rates fin, fhi,fh3, P12 for the i th 
node are chosen to be larger than respective parameters in the f h node (i > j, i, j = 1 , 2, 3, 4). The reproduction number increases 
with a) 2 v livestock movement rates from node j to node i, increasing in each run and decreases when the movement rates from 
node i to node j increases as is shown in Figure [2(b)| In scenario 2, where 0*2; > ^2j, the reproduction number decreases when 
the movement rates from node j to nodes i increase and increases when the movement rates from node i to nodes j increase as 
is in shown in Figure [3(a)| and Figure [3(5)1 respectively. When livestock recovery rates 72/ > J2j in scenario 3, the reproduction 
number decreases with the increase of movement rates from node /' to node i and increases with the increase of movement rates 
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(a) The reproduction number and its lower bound with (b) The reproduction number and its upper bound with 
homogeneous parameters. homogeneous parameters. 

Figure 1: The reproduction number and its upper and lower bound for 100 nodes with homogeneous parameters except movement rates. 



No. 


parameter 


movement rates 


Ro 


1 


Pm > Puj, fan > Puj, Pm > Pi3j, Pm > p-iij 


a>ji increases 


increases 






a>jj increases 


decreases 


2 


U2i > C$2 j 


tJjj increases 


decreases 






a>ij increases 


increases 


3 


J2i > 72 j 


ojji increases 


decreases 






oiy increases 


increases 


4 


m > H2j 


ojji increases 


decreases 






a>ij increases 


increases 



Table 2: Different scenarios for numerical simulation on four node networks. The superscripts j = 1, 2, 3, 4 and i > j. 



from i to j as shown in Figure [4(a)| and Figure [4(b)! respectively. Similarity, when livestock mortality rates /i2i > H2j in scenario 
4, the reproduction number decreases with the increase of moment rates from node j to node i and increases with the increase 
of movement rates from i to node j as shown in Figure |5(a)| and Figure |5(b)| respectively. In all scenarios, Rq can fluctuate 
from below 1 to above 1 or from above 1 to below 1 by tuning the parameters. Figure [3] Figure [4] and Figure [5] have shown a 
similar trend which is opposite to that of Figure [2] The simulation results of the above four secnarios have shown that livestock 
movement rates affect the reproduction number when there is spatial heterogeneity. 




increase when /?2, > Pij Ro increases. increase when fix, > fiij, Ro decreases. 

Figure 2: The reproduction number for four nodes with different contact rates. 
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1.005 
1.004 
1.003 
1.002 
1.001 
1 

0.999 
0.998 







20 40 60 

Simulation index 




20 40 60 80 

Simulation index 



(a) As the movement rates from node j to node i (or-j) (b) As the movement rates from node £ to node j (ojfj) 



increase when d2i > chj, Ro decreases. 



increase when dn > d2j, Ro increases. 



Figure 3: The reproduction number for four nodes with different livestock death rates. 




Figure 4: The reproduction number for four nodes with different livestock recovery rates. 




Figure 5: The reproduction number for four nodes with different livestock mortality rates. 



4. Results and discussions 

We have proposed a general multi-species model with both vertical and horizontal transmission. The formula of the value 
of the reproduction number has been derived for h species, including g species vertically transmitting the disease and h — g 
species transmitting the disease both vertically and horizontally. The expression of Ro is a function of vertical and horizontal 
transmission as is shown in Equation d38l >. 
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For Rift Valley fever meta-population model, the expression of Rq and its upper and lower bound have been derived. The 
reproduction number can be estimated by lower and upper bound if its exact value is hard to obtain due to the size and the number 
of species involved in a network. The difference between the upper bound and the lower bound of the reproduction number is the 
maximum transovarial transmission rate qu across all nodes as is shown in Equation (f41~b . Therefore, the tightness of the bounds 
is determined by the largest transovarial transmission rate on the network if the reproduction number for horizontal transmission 
R^ is computed by Equation d47b . If the reproduction number for horizontal transmission is to be estimated, the tightness of the 
bounds for Rq depend on disease transmission parameters and movement rates between different locations according to Equation 
(SB and @5). 

If all nodes are spatially homogeneous with respect to disease parameters, the reproduction number is independent of the 
number of nodes, the size of each node, and the movement rates between different node shown by simulation results. The bounds 
of Ro is proved to be independent of movement rate in Theorem (0). Moreover, the expressions of the lower and upper bound 
of Rq are the same as those presented for homogeneous populations without considering movement presented by 1)3 ill . The 
reproduction number for horizontal transmission, R^, shown in Equation (T59t is proved to be only related to disease parameters 
for homogeneous networks. 

We have studied the conditions when movement rates are playing a role on Rq. We have observed consistent trends of 
the reproduction number when we change some movements rates between nodes with smaller values to larger values of some 
parameter or in the opposite direction. Rq can fluctuate above or below 1 . The observation helps us to design efficient mitigation 
strategies by executing movement ban between some nodes and in some directions. Therefore, livestock movement plays an 
important role on the dynamics of epidemic spread. The simulation results have shown that movements among different nodes 
only affect the reproduction number when the network is spatially heterogeneous. Livestock death rate, incubation rate, mortality 
rate, recovery rate, and contact rate can be different for livestock in different nodes due to weather conditions, public health 
facilities, environment, or the type of nodes, such as feedlots with higher death rate and livestock premises with lower death 
rate. Different weather conditions such as rainfall affects the birth rate of mosquitoes, temperature affects mosquito incubation 
rate. Even if weather conditions are homogeneous across all nodes, different genera or species of mosquitoes may have different 
incubation rates, contact rates, or birth rates. Whatever the heterogeneity between nodes is, the same mathematical model, the 
expression of the exact reproduction number and its bounds still can be applied. The expressions of the reproduction number 
and its bounds presented in this paper can be easily applied to many other disease models besides the Rift Valley fever model. 
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Theorem 5. Let A k (k — 1,2,3, • • • ,m) be all n X n invertible matrices such that A k are positive or irreducible nonnegative 
matrices. Denote the (i,j) entry of A k by atij. Let a k — min,(2i fly) > 0, a k — max;(Zi a ij)> men I~KLi ~n ^ pCil^i^ 1 ) ^5 



Appendix 




Proof Clearly, 



a\C < CA k < agC, 



where 



C= [1 1 1 1], 



Since A~ l ^ 0, 



a\C ^ CA k => a L k CA k l sC CA^ 1 => CA A T [ sC \c => CA k l A k \ ^ —CAl\ < 




1 



C 




By repeating this process, we have 




Recall, for any nxn positive or irreducible nonnegative matrix A, the following property is satisfied 114 111 . 



n 



n 




Because the entries of C Yl 




is the sum of each column of matrix Yl' 




by the above property, 
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Theorem 6. Mi and Xz in Vh are invertible. Moreover, Matrix M 2 1 and X 2 1 are nonnegative matrices. Matrices M2 and X2 
are matrices in <\39\ and <\40\ , respectively. 

Proof Note that M2 is a diagonal dominant matrix of its column entries. The proof that M2 and X2 are invertible is omitted 
because it is similar to the proof for proving 'W is invertible in TheoremQ] Next, we prove M 2 l is nonnegative. Matrix M2 can 
be rewritten as follows. 



Mi 



d 2 \ + £21 + 2" =2 0) 
- ,2 



"12 
,2 

' J l» 



dz2 + £22 + H" = ij # 2 ( 
- ,2 



'2; 



,J 2n 



-or, 
-or. 



din + B2n + Z"=l 0) 2 nj 



(68) 



dz\ + £21 + Yl)=2 ^\ 








• • • d 2n + B2„ + Z"=l <^ 2 „j 

G-H 






W 21 ' 




W ?2 
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■ Kz 


< 




■ 



Therefore, 



and 



Moreover, 



! fv 21 , - , V" / .2 
— +E!+i i= 2 M lj 






1 



ff 2 z ^J=l «„ 





"21 

■ / : A ':i , „ , v" , ,2 
n^- +e 2+2, J= 2 w l; 





^2 ^.7=1 nj 



< J](G- l H) u < 1 

;'=i 








■^21 , „ , v» , ,2 



Therefore, G H is convergent. Obviously, G 0, and G H ^ 0. Equivalently, M2 is a M-matrix and M 2 ^ [42]. 
Similarly, X2 is a M-matrix and Xz 1 ^ 0. 
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